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ABSTRACT: Hawking’s model of black hole evaporation is not unitary and leads to 
a mixed density matrix for the emitted radiation, while the Page model describes 
a unitary evaporation process in which the density matrix evolves from an almost 
thermal state to a pure state. We compare a recently proposed model of semiclas- 
sical black hole evaporation to the two established models. In particular, we study 
the density matrix of the outgoing radiation and determine how the magnitude of 
the off-diagonal corrections differs for the three frameworks. For Hawking’s model, 
we End power-law corrections to the two-point functions that induce exponentially 
suppressed corrections to the off-diagonal elements of the full density matrix. This 
verifies that the Hawking result is correct to all orders in perturbation theory and also 
allows one to express the full density matrix in terms of the single-particle density 
matrix. We then consider the semiclassical theory for which the corrections, being 
non-perturbative from an effective field-theory perspective, are much less suppressed 
and grow monotonically in time. In this case, the Renyi entropy for the outgoing 
radiation is shown to grow linearly at early times; but this growth slows down and 
the entropy eventually starts to decrease at the Page time. In addition to compar¬ 
ing models, we emphasize the distinction between the state of the radiation emitted 
from a black hole, which is highly quantum, and that of the radiation emitted from 
a typical classical black body at the same temperature. 
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1 Introduction 

There has been a recent spike in activity on understanding the implications of black 
hole (BH) evaporation [1, 2]. This can be attributed, in large part, to the contro¬ 
versial proposal that a unitary evaporation process comes with the cost of a “fire¬ 
wall” [3] or “energetic curtain” [4]— these being colorful euphemisms for an apparent 
tension between general relativity and the unitarity of quantum field theory in a BH 
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setting [5-8]. While the debate rages on, a consensual mechanism for information 
release is still lacking. 

In Hawking’s model, the process of BH evaporation is not unitary [2], Hawk¬ 
ing argued that the correlation functions of the emitted radiation are diagonal in 
mode-occupancy number, frequency and emission time and, from this, deduced that 
the density matrix for the radiation is diagonal in the same quantities. In fact, 
to good accuracy, the final state of the emitted radiation is thermal, and so it is 
similar to a maximally mixed state. Then the evaporation cannot be the result of 
unitary evolution from a nearly pure state, which lead Hawking to state (a statement 
which he more recently retracted [9]) that the process of gravitational collapse is not 
compatible with the standard principles of quantum mechanics. 

The established benchmark model for describing unitary BH evaporation has 
been the Page model [10]. In this model, the BH and the emitted radiation are 
assumed to be in a pure state in some large Hilbert space, which is partitioned 
into an “in” part representing the BH and an “out” part representing the outgoing 
radiation. Because the combined state is pure, there must be a special basis in which 
the full density matrix has only a single eigenvalue. An external observer measures 
the reduced “out” density matrix in some random basis and treats the transformation 
matrix U as a random matrix with prescribed statistics [10, 11]. (See also [12] for 
a recent review.) When the “out” system is smaller than the “in” system, it looks 
to an external observer as if it were thermal. Conversely, when the “out” system 
becomes the larger one, the deviations from a thermal state grow and indicate that 
the total state is indeed pure. The critical time in which the midpoint of evaporation 
is reached is normally called the Page time. Given that the BH is initially in a pure 
state and that the evolution is unitary, the Page model can plausibly be viewed as 
setting the minimal rate for purification of the emitted radiation. 

Two of the current authors (RB and AJMM) have recently proposed a semi- 
classical model of BH evaporation [13, 14], with the premise of repeating Hawking’s 
seminal calculations [1, 2]— which assume a classical background metric—so as to 
include a fluctuating BH geometry. The incipient BH is endowed with a quantum 
wavefunction [15, 16], leading to expectation values in place of fixed classical param¬ 
eters. From this perspective, Hawking’s model pertains to the limit of an infinitely 
massive BH with a fixed size, whereas the semiclassical reformulation treats the BH 
mass as finite with a continually decreasing size due to classical back-reaction effects. 
In the new model, the correlation functions of the emitted radiation are no longer 
diagonal and the evaporation process becomes unitary, even though the thermal-like 
emission spectrum is maintained. This model of unitary BH evaporation obeys the 
general constraints about the rate of purification. In particular, it has been shown 
in [17], relying on the results of the current paper, that the Renyi entropy of the 
radiation decreases monotonically from its peak value at the Page time and does so 
at a faster rate than the Page model predicts. 
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From the point of view of an effective field theory in a fixed curved-space back¬ 
ground, the correlation functions of the semiclassical model are modified in a non- 
perturbative way [16]. Contrary to expectations, these modifications could become 
significant in some situations. It has been shown that, when the leading result in 
the effective field theory either vanishes or diverges, then the non-perturbative cor¬ 
rections are particularly relevant [18]. In the first case, we expect a correction of 
the order 1/Sbh (the inverse of the BH entropy) to become the leading result and, 
in the second, the divergence should be replaced by a large but finite Sbh ■ These 
small non-perturbative quantum effects could be coherent, and so their amplitudes 
could add up and grow with time; even growing to the point where they become 
comparable to the classical outcome [13]. 

The previous studies of the new model have utilized the single-particle density 
matrix (the two-point function) of the radiation as the primary tool 1 rather than 
the full multi-particle density matrix. The single-particle density matrix is easier 
to calculate and, for the radiation emitted by a BH, sufficient to completely deter¬ 
mine the full density matrix. One of the primary goals of the current paper is to 
make explicit the relationship between the single-particle and multi-particle density 
matrices. The results are presented in a way that can be applied to the Hawking 
model, the proposed semiclassical model and even to more general modifications of 
Hawking’s framework. As an upshot, we are able to find a closed expression for the 
Renyi entropy of the full density matrix that is expressed directly in terms of the 
single-particle density matrix. 

It has been argued for BHs in Anti-de Sitter (AdS) space—at first by Malda- 
cena [19] with subsequent elaborations by many others — that non-perturbative 
contributions from other geometries which do not possess a horizon could be rele¬ 
vant to the state of the Hawking radiation. Maldacena’s specific example was the 
contribution of thermal AdS to the matter correlations for large BHs. In this case, 
the contribution to the density matrix of the radiation is expected to be exponen¬ 
tially suppressed ~ e~ SsH . However, because the number of off-diagonal elements 
is exponentially large ~ e +SBH , one could not be sure as to the relevance of such 
a contribution. It could then be argued that, without a better control on the non- 
perturbative contributions from other geometries (which are usually not available), 
it becomes difficult to trust the Hawking result. 

We have two comments on this issue. First, according to our analysis, one 
requires much larger corrections ~ e ~ SBH / 2 to challenge the reliability of the leading- 
order outcomes (See, for example, [20]). Second, the quantum fluctuations of the 
background BH geometry induce off-diagonal effects in frequency space that are only 
power suppressed ~ 1/ \/Sbh and sufficient by themselves to restore unitarity. These 

^^This matrix was called “the density matrix” by abuse of language or “the radiation matrix” in 
the previous studies. 
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statements will be made more precise in due course. 

The rest of the paper is organized as follows: In the next section, we review 
the basic framework of the Hawking model of BH evaporation and recall its main 
results. Then, in Section 3, we formulate a precise relationship between the multi¬ 
particle density matrix and the single-particle density matrix. The preceding analysis 
is applied to the Hawking model in Section 4, with the aim of regularizing formal 
divergences that arise in the eternal-BH framework. Section 4 also includes a dis¬ 
cussion on how the state of the BH radiation differs from that of thermal radiation 
emitted by a “normal” black body. The formal analysis up to this point sets the 
stage for Section 5, where we discuss similarities and differences between the three 
models of BH evaporation: the Hawking model, the Page model and the semiclassical 
model. There is then a brief conclusion in Section 6, followed by four appendices 
that fill in some technical gaps in the main text. 


2 The Hawking model of black hole evaporation 

Let us start here by recalling the original Hawking description of BH evaporation, 
which completely dismisses the back-reaction and time-dependence effects. This 
paper uses the notations of Hawking in [2], In particular, we consider the initial 
vacuum state |0_) and a final vacuum state |0 + ). The initial vacuum contains no 
particles with respect to the creation and annihilation operators {a+, a;} as defined 
at past null infinity, 

Oi|0_) = 0. (2.1) 

Following Hawking, the Hilbert space at future infinity is a direct product of 
the Hilbert space of the particles falling into the BH, H- m , and of those going out, 
'Hont, so that % = "Hm <S> Tfout- The ingoing particles are created and annihilated 
by the operators {c+, c,}, whereas the outgoing ones are created and annihilated by 
{bf, bi}. Both are related to { af , a;} via Bogolyubov transformations. The final 
vacuum state can be represented as |0+) = 10i n ) 10 ou t) and contains neither outgoing 
nor ingoing particles, 

6i|0+) = 0, c i |0+) = 0. (2.2) 

The two vacuum states |0_) and |0+) do not coincide with one another due to 
the particle creation by the BH. In other words, an initially empty state will appear 
to have a non-zero occupation number in the final state, meaning that 6j|0_) ^ 0. 
The initial vacuum can be expressed as a linear combination of the Fock states at 
future infinity as 

|0_) = A a6 |a out )|6 in ) , (2.3) 
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where |a out ) and |6 in ) are the basis vectors of the Fock spaces which are spanned by 
the operators {bf, bi} and {cf, Ci}, respectively. As a consequence, 


k„ t > = n —j=i (k)” is l°»"*> • IU = n -r=, (4)”“ | 0 ,„). (2.4) 

j \Z n ja- V Vrikb'- 

The object of interest is the density matrix of the initial vacuum state, 


p vac = |0_)(0_| , (2.5) 

which allows one to calculate the vacuum expectation values of observables according 
to the standard rule (0_|C?|0_) = tr p va , c O. 

In particular, we will be interested in the expectation values of the observables 
at future infinity. These will be composed only from the outgoing creation and 
annihilation operators {bf, bi} and so can be written as O — O orit <E> X in , where 
X in denotes the identity operator in the Hilbert space of the ingoing states. The 
expectation value of this operator then takes the form 

<0_|O|0_) = <0_|a out |0_) = p out ac O out ca = tr out p out a out . (2.6) 

Here, p out ac are the matrix elements of the reduced density matrix of the Hilbert 
space "Hout for the outgoing radiation and are determined by 

P ° ut = tr in p vac = (6 in |0_)(0_|& in ) = A ab A cb |a out )(c out | = p out ac |a out )(c out | , (2.7) 


whereas 

O ont ca = (cont\d out \a out ) . ( 2 . 8 ) 

The elements of the density matrix for the outgoing radiation can be found by 
calculating different moments of the operators {bf, bi}. In particular, the number 
operator of the jth outgoing mode is 

(rij) = (0 _|bf bj |0_) = P™*™ ■ ( 2 - 9 ) 


The famous result of Hawking is that these elements are the same as for thermal 
radiation, 


P 


out 


fin ja nj c P{ n ja) ) 
j 


Pin) 


(1 - e- w/T )(e- a;/T r) n 
[1 - (1 - Y)e~ UJ / T ] n+1 ’ 


( 2 . 10 ) 


where T is the grey-body factor. The density matrix (2.10) is diagonal in both j a 
labeling the mode frequencies and Uj a , the occupation number of each mode. The 
index a (also c) runs through the set of basis vectors of the Fock space that is built 
by the operators bf. Hence, the dimension of the density matrix is formally infinite. 
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The result of Hawking indicates that the density matrix for the outgoing radi¬ 
ation is thermal and thus is similar to that of a maximally mixed state. This is in 
contradiction to unitary evolution and at the core of the information paradox since 
the collapsing matter was in a pure state and evolved into a mixed state. Nonethe¬ 
less, we will eventually show how the presence of off-diagonal elements in the density 
matrix p out can lead to the purification of the outgoing radiation in the course of 
the BH evaporation. 


3 The multi-particle density matrix in terms of the single¬ 
particle density matrix 

The goal of this section is to expose the connection of the full density matrix of the 
outgoing radiation p out , as defined in (2.7), to the widely used single-particle density 
matrix 

p‘l = (0_|6+6‘|0_> . (3,1) 

Let us begin with the well-known Bogolyubov transformation between the cre¬ 
ation and annihilation operators at future infinity { bf, b *} and at past infinity {a^~, a 1 }. 
This is expressible as 2 




(3,2) 


In terms of the Bogolyubov coefficients, we can express p® ■ as follows: 


(3.3) 


As shown in Appendix A, there is a closed expression for the density matrix p° ut 
in terms of the V s, 

-out = I e -b+nrv ; (3.4) 

Z 

where Z is the normalization factor and is a c-number Hermitian matrix. The 
latter is related to the single-particle density matrix ph as 


1 n 3 1 

p = -- , or e = 1 + - . 


& 1 — 1 p 

From now on, the superscript “out” on p will be omitted for brevity, 


?=r nt 


(3.5) 


(3,6) 


2 Here, upper and lower indices have been introduced so as to make the matrix multiplication 
explicit. The summation symbol will typically be omitted in what follows. 
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The full density matrix p is thus completely defined by ph = (3 lk (3jk, with the 
normalization factor Z expressed as 


Z 


det 


1 

1 — e~ n 


det [1 + p\ . 


(3.7) 


One can define the entropy and related quantities directly in terms of the matrix 
p. In particular, an explicit expression for the Renyi entropy, H 2 = -ln^CJ, is 

(tr p) 

readily obtained by evaluating 


tr [p 2 ] =^ tr 

= det 

from which it follows that 


-2 b+iib 


1 — e 


-n 


1 + e 


-n 


= det [l — e n ] det 
= det 


1 — e 


-20 


1 + 2p 

H 2 = tr [In (1 + 2 p)\ . 


(3.8) 


(3.9) 


From this expression, it would appear that the Renyi entropy does not vanish 
unless p = 0, when no particles have been emitted. However, the entropy for a state 
with N radiated particles can reach as low as H 2 = In (2 N) <C Sbh (0), which is below 
the scale of validity of the approximations that were used in deriving Eq. (3.9). 

The expression for the von Neumann entropy is more cumbersome but we include 
it for completeness: 


S 


tr [pin p] = tr In [1 + p] + det 




(3.10) 


In this paper, we will therefore concentrate on the Renyi entropy H 2 . It is much 
easier to calculate and provides an accurate measure of the entanglement in the 
state of the BH radiation. 3 


4 The state of the emitted particles in the Hawking model 

Having a formal expression for the density matrix at hand, we can now evaluate the 
Renyi entropy and the particle number in the state of the outgoing radiation in the 
Hawking framework. However, as made clear below, these quantities are formally 
proportional to an infinite sum over all the possible modes in the Fock space of the 
outgoing radiation. This divergence comes about from implicitly assuming an eternal 
BH that radiates for an infinite period of time and, hence, emits an infinite amount 
of particles. 

3 For Gaussian states, such as the vacuum state of a non-interacting theory, H 2 respects the 
strong-subadditivity condition in the same way the von Neumann entropy does [21]. 
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The physical situation, however, must be different: The number of particles 
emitted by the BH during a finite time interval is finite and depends on the initial 
mass of the BH. The goal of this section is to present a way of counting the modes 
emitted by the BH which reach an observer at future infinity during some finite time 
interval. 4 By implementing this method, we find that the physical quantities indeed 
become finite. For simplicity, the analysis is initially specialized to the case of the 
Hawking model. 

In the current case, the matrix Qh = ^<5*■ ■ is diagonal, and Eqs. (3.7) and (3.9) 
for the partition function and the Renyi entropy then read 


^ln(l-e-^ T ) , 

(4.1) 


(4.2) 


while the total particle number is given by 

N = e u>i/T _ i ■ ( 4 - 3 ) 

i 

The sums in the above expressions run over the full set of the frequency modes 
uji of a quantum field in the BH background. Hence, the dimension of the set of the 
possible frequencies is infinite, leading to an infinite entropy and particle number. 


4.1 A reminder: A scalar field in a box 


To better understand the correct way of counting the modes so as to render physical 
quantities finite, let us first recall a simple case: A massless quantum field in a box. 

A massless quantum scalar field 0(£, x) in a box of size L with reflecting boundary 
conditions (he., </>(t, 0) = 0(£, L) = 0) can be expanded in terms of the Fourier modes 


u^(t,x) 


1 1 

-sin 

(2ttL) 3 / 2 (2cu)V 2 



sin 



sin 



(4.4) 


with frequencies u = ^ \Jk\ + /c 2 + /c 2 , where k x , k y and k z can be any integer. 

For a large enough box (L 1/T), the sum over the modes can be approximated 
by an integral over a continuous momentum p, 


£ = W 

L. 7 „ ’’ 


d 3 p 
( 2^)3 ’ 


(4.5) 


with the corresponding frequency given by ui(p) = \p\. In this way, the quantities in 
Eqs. (4.1)-(4.3) are proportional to the volume of the box V = L 3 ; that is, they are 

4 For an earlier, detailed calculation of the physical emission rates that employs a different tech¬ 
nique, see Ref. [22]. 








extensive. Moreover, since the only other dimensional parameter is the temperature 
T, any extensive, dimensionless quantity is given by V T 3 times some numerical 
constant: 

1 n 2 

-lnZ = VT 3 - — dxx 2 In (l — e~ x ) = ——V T 3 , (4.6) 

27 t 2 J 0 90 

ft = K T* • ^ (i±^) = ^ TV (4.7) 

N = VT 3 -^r dx ^i = ^ VT3 ’ < 48 > 

where x = ui/T. 

In general, due to the diagonal form of the thermal density matrix, any quantity 
that is given by a trace over the Hilbert space % out is proportional to a product over 
the available frequencies (as in Eq. (A.7)). A logarithm of such a quantity is then 
proportional to a sum over the frequencies and is extensive by virtue of Eq. (4.5). 
Imposing some physical boundary conditions on the scalar field results in a restriction 
on the allowed frequency modes—in this case, the condition that the wave numbers 
have to be integer multiples of 27 t/L. The formally infinite quantities (4.1) - (4.3) 
then become finite, as should be true for any physical system. 

4.2 Emission of localized wave packets 

We have just seen that, for a held in a box, the frequencies are quantized in units 
of the inverse size of the box and the total number of modes is proportional to the 
volume of the box. In order to study how a BH radiates into an infinite space, it is 
more convenient to use localized wave packets with finite normalization instead of 
the Fourier modes (as in [1]). 

Let us begin here with a scalar held that is propagating in the background of 
a BH. It can be expanded in terms of the complete set of solutions of the wave 
equation having only positive frequencies. At past infinity, the expansion contains 
only ingoing modes and takes the form 

0 = 5^ + M + ) ’ ( 4 - 9 ) 

i 

where {a) 1 ", a*} are the creation and annihilation operators defined at past infinity and 
the sum runs over a discrete set of modes {_/)} with finite normalization. At future 
infinity, the expansion of the scalar held in terms of positive-frequency solutions 
contains both ingoing and outgoing modes, 

0 din T bout j bout ^ ^ ( Pi^i pibj ) , (4.10) 

i 

where {bf, bi} are the creation and annihilation operators of the outgoing modes 
{pi} at future infinity. 
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In the continuum normalization with frequencies oj (rather than the discrete wave 
packets i), the ingoing and outgoing Fourier modes for the spherically symmetric 
solutions are expressible as 


fuim{v,r,9,<f>) = F ulm (r)Y l m (9,cf ) )e iu}V , (4.11) 

Pu,im{u,r,9,<j>) = P ul i m {r)Y l m {9A)e^ u , (4.12) 

where l and m are the angular-momentum numbers, and v and u are the advanced 
and retarded tortoise coordinates (respectively). 

A complete set of localized wave packets can then be defined as [1] 

rU+ i)e 

fjnlm(v , T, 9, 0) = £~ 1/2 / e -2mnu>/e 0) <2CU , 

Jje 

rU+ 1 ) £ 

Pjnim(u, r, 9, 0) = e _1/2 / e~ 2irinu}/£ Puimiu, r, 9 , 0) dcu . 

Jje 

In place of the continuous label u, the wave packets have been labeled by two integer 
indices j > 0 and n. The wave packet with index j contains waves with frequencies 
that are localized in the range j £ < co < (j +1) £. Index n labels the ray along which 
the packet is propagating: The wave packet fj n is peaked around the advanced time 
v = 2nn£~ 1 and the wave packet pj n , around the retarded time u = 2i me -1 , both 
having a width of 27r£~ 1 . 

Our particular interest is the wave packets which were emitted by the collapsing 
body and then detected at some fixed distance away from the BH during a given 
period of time t 0 < t < t 0 + A t. Such wave packets will also be localized in the 
corresponding range of retarded time Uq < u < Uq + At. For a long enough time 
interval, one can always choose the width e -1 so as to ensure the detection of many 
wave packets (e -1 At), as well as maintain a fine-enough frequency resolution 
(e T) for each wave packet to be treated as a monochromatic mode of fixed 
frequency. In such a case, the summation over the wave-packet position n in the 
interval An = At/(27re _1 ) can be approximated as 


(4.13) 

(4.14) 


no+An 

E 

n=riQ 


2vr 


At, 


(4.15) 


whereas the summation over the discrete set of frequencies j = u /e can be approxi¬ 
mated by the integral 

°° r du 


E 

j =o 


(4.16) 


As one can now see, in this approximation, the total number of modes which 
can be detected during the time interval At, does not depend on the choice of the 
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parameter e and is proportional to At times the integral over mode frequencies, 


oo no+An 

E = E E 

i j =0 n=no 


At 


r°° du 
Jo 2vr ■ 


(4.17) 


The total number of particles detected during the time interval At in this case 
is given by 


oo no+An oo l -p 

N = H E ZE+t 

j =0 n=no 1=0 m=—l 


oo l 

4, EE 


0 m=—l 



dtd r ujlm 
2ir e u l T — 1 ’ 


(4.18) 


where T^ m , the so-called grey-body factor, is determined by the properties of the 
modes near the horizon. 

The emission of the modes with high multipoles is highly suppressed, T,^ m <C 
T w0 o = r(tu) for l^m > 0, and the sum (4.18) can be well approximated with the l = 
m = 0 term only (see, e.g., [22]). In what follows, we will consider only the l = m = 0 
modes and omit the angular-momentum labels. Moreover, in the Schwarzschild case, 
there is only one available dimensional parameter, which can be chosen as the BH 
temperature T. Therefore, the grey-body factor T(o;) depends on the frequency only 
through the ratio x = uj/T, and the number of emitted particles simplifies as follows: 


N 



dx T(x) 
2tt e x — 1 


r\j 


At AT 3 



dx T(x) 
2ir e x — 1 ’ 


(4.19) 


where A ~ T 2 is the area of the BH horizon. 

Thus, the rate of BH emission coincides with the thermal emission rate of a body 
with area A. In reality, due to the non-trivial frequency dependence of the grey-body 
factors, the spectrum for BH radiation is quite different from the thermal spectrum 
of a body in an empty box. Some details about the grey-body factors are provided 
in Appendix B. 

And so, in the absence of other dimensional parameters, all the extensive quanti¬ 
ties describing the emitted radiation are proportional to the product At T or, equiv¬ 
alently, to the total number of emitted particles (4.19). For example, the Renyi 
entropy H 2 is given by 


H 2 


At 




awm 

e oj/T _ 1 J 





oc N . 


(4.20) 


Hence, by properly accounting for the finite duration of detection, one finds that 
the number of particles emitted by the BH, as well as the entropy of the radiation, 
becomes finite as opposed to the idealized case of infinite emission time. 
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4.3 Single-particle density matrix in the wave-packet basis 

In order to highlight the connection between the localized wave-packet description 
and the semiclassical model as presented in [13, 14], we shall construct the single¬ 
particle density matrix p l - in the finite wave-packet basis. 

Using the definition of the wave packets (4.14), one can express the matrix 
elements of p in terms of the same Fourier-mode basis, 


in 

P j'n' 


= £ 


-1 




dcue 2 ™ /£ 


hi'+ik 


dw' e- M “' /£ p(uj,uj') , 


(4.21) 






where 

/•OO 

p(u, u/) = dti PuuPu'u ■ (4.22) 

Jo 

Each wave packet is labeled by the characteristic frequency uij — j £ and the 
detection time t n = 27r£ _1 n. The choice of £ parametrizes a trade-off between the 
frequency and the position resolutions. The total number of modes W, within a given 
range of frequencies Au; and an interval of time At, does not depend on £ and is just 
given by W = A j ■ An = At Au>/ (2tt) . In the case of a thermal distribution, the 
frequency range Au; is effectively set by the temperature. Hence, the total number of 
modes with substantial occupation numbers is of the order W = T At. The matrix 
P >n j'n' Eq. (4.21) is, therefore, effectively a W x W matrix. 

An important feature of BH radiation is that the emission rate T is the same as 
the temperature T. For a general radiating body, this is not true; even if T oc T, the 
proportionality constant need not be close to unity. This has important consequences, 
which we now elaborate on. 

From the above discussion on wave packets, we have learned that the total num¬ 
ber of emitted modes is given by W = AuAt, where Au; is the range of emitted 
frequencies and At is the detection time. For the case of a nearly thermal emitter, 
the frequency range is given by the temperature Au; ~ T, so that W = TAt. On the 
other hand, the total number of emitted particles during this same time is given by 
N = TAt, where T is the emission rate. For a BH, with the emission rate T ~ T, the 
number of emitted particles N is then of the same order as the total number of oc¬ 
cupied frequency modes W; i.e., N ~ W. This fact has two important consequences 
for BH radiation: 

First, the average occupation number for the modes of radiation is approximately 
of order unity, as N/W ~ F /T ~ 1; meaning that, on average, each mode of radiation 
is occupied by only a few particles. Second, since W ~ N, the single-particle density 
matrix p is effectively of size N x N. Then, as the average occupation numbers 
correspond to the diagonal elements of the matrix p, these elements go as (9 n - n ~ 1. 
Hence, for BH radiation, the two-point function p in the wave packet basis has the 
form of an N x N identity matrix. 
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The fact that, on average, each mode is occupied by only a few particles also 
implies that the radiation field in a BH background cannot be treated as semiclassical. 
This is not the normal state of affairs. For any other thermal emitter, the radiation 
rate is proportional to its area, T ~ A ■ T 3 , and the area for classical emitters is 
much larger than the inverse temperature squared (i.e., the square of the typical 
wavelength), so that T/T ~ AT 2 1. Consequently, N/W ~ T/T » 1 and the 
radiation field is highly classical. 

In order to make these ideas more precise, let us calculate tr p and tr p 2 for the 
Hawking model [1, 2], For a large BH, 

P+, u') = 5{uj - u') p(u) = S(u - u') x . (4.23) 

With this expression, in the wave-packet basis, 

rU+ 1 ) s 

fA n j'n> = £ ~ 15 ]’ / duj e2 ™ (n ~ n } w/£ P M • (4.24) 

As one can observe, p jn 3 , n , depends only on the difference in detection times 
t n — t, n i. Because of the finite time and frequency resolution, it has non-zero off- 
diagonal elements with n ^ n', which are concentrated in a narrow strip \n—n'\ < e/T 
corresponding to \t n — t n '\ < T~ l . The number of non-zero elements of p is, therefore, 
proportional to its dimension W and, consequently, of the order of the number of 
emitted particles N. 

The traces of p and p 2 can be explicitly calculated: 


oo no+An 
j =0 n=no+l 


in 

P in = £ 


no+An roo roo ^ 

£ / du p(u) = At — p{u) ~ At T 

n=n 0 + iV 0 ^0 71 


W ~ N , 
(4.25) 


oo no+An 

i 2 \ ^ \ ^ jn j'n' 

tr P =£ £ P j'n'r jn 

j= 0 n=no+l 
j '=0 n'=no+l 

no+An r(j+ l)e 


= £ 


-2 


£ E 

j =0 n=no+l 3 £ 
n'=no~\-l 


du / du' p(u) p(u') 




= £ 


-2 


£ 


'■0'+ 1 ) E r(j +!) £ 


du 


du 


At 


3=0 J 3 £ 

f°° du 




, ( sin [An n (u — u') / e\ 
\ sin [ 7 r + —u/)/e] 


p(u) P +') 


'0 


2n 


p{u') 2 ~ At T ~ N 


(4.26) 


In order to obtain the last line, we have used that (sin(Amrx)/sin+x)) 2 ~ A nS(x) 
for large An. 
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The number of the non-zero eigenvalues, which is just the number of diagonal 
elements in the Hawking model, can be then estimated by the participation ratio 
(i.e., the inverse of the purity; see, e.g., [23]), 


PR = 


(tr[p]) 2 

tr [ p 2 ] 


(4.27) 


The result grows linearly with the number of emitted particles and also provides an 
estimate of the Renyi entropy (4.20). 

For a more general case including off-diagonal elements, one can still expect 
tr p ~ N, as this result only depends on the dimensionality of the matrix. On the 
other hand, tr p 2 can be much different, but its value can be readily estimated if the 
number of non-zero eigenvalues is known or vice versa. 


5 Off-diagonal elements of the density matrix 

With the previous framework in hand, we will now consider the structure of the 
off-diagonal elements in both the single- and multi-particle density matrices. This 
will be carried out for the three distinct models of an evaporating BH: (1) Hawking’s 
model [1], (2) the Page model [10] and (3) the semiclassical model [13, 14], 5 

In this section, all dimensional quantities are expressed in units of the BH 
temperature, N = N(t ) is the total number of emitted particles at time t and 
Sbh = Sbh{N ) = S bh (0) — N{t) is the BH entropy at “time” N. This relation 
is modified, semiclassically, by fluctuations in the emission time of the Hawking 
modes. However, the fluctuations were shown to be small [13], being proportional to 
inverse powers of Sbh- 

One of the main issues to be clarified is to what extent the off-diagonal elements 
in terms of mode-occupation number of the full density matrix are suppressed. If 
they are indeed highly suppressed, then the relationship between the density matrix 
p and the single-particle density matrix p is reliable. 

5.1 The Hawking model 

In Hawking’s model, the outgoing radiation is described to a very good approxi¬ 
mation by the diagonal density matrix given in Eq. (2.10). What we would like to 
address here is how accurate an approximation this is by quantifying the off-diagonal 
contribution. 6 We shall, however, not discuss any back-reaction nor evolution effects 

5 For the Page model, the single-particle density matrix cannot be calculated without additional 
information because this framework is based upon choosing an arbitrary random basis in the com¬ 
bined radiation and BH Hilbert space. How this basis is related to the Fock space of the radiation 
is left unspecified. 

6 By off-diagonal contribution, we mean that beyond the broadening effects of the wave-packet 
treatment in the last section. 
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to the Hawking model in this subsection. These effects were originally disregarded 
by Hawking himself and introducing them would mean deviating from the original 
model in a significant way. We do account for the evolution effects in the semiclassical 
model that we discuss in the subsection 5.3. 

The Hawking density matrix p# is only approximately diagonal in two obvi¬ 
ous ways: Elements between states differing in total occupation number are only 
approximately vanishing, as are those between states with some modes differing in 
frequency or in the time of emission. We would like to estimate the magnitude of 
the off-diagonal elements of the density matrix. As far as we know, this issue has 
not been addressed directly in the literature, although Hawking did expect such 
corrections to be exponentially small [1], 

What we find below is an overall suppression for off-diagonal elements in mode- 
occupation number that is stronger than e~ SBH . This is significant because it implies 
that the results of Hawking are not modified in a meaningful way by off-diagonal 
corrections.' We also find that the off-diagonal elements of the single-particle density 
matrix p l 3 are smaller than 1/Sbh, leading to corrections to physical quantities, 
such as the entropy, that are suppressed by at least a factor of 1 /S 2 BH . These are 
negligible at all times, as originally anticipated by Hawking. Corrections to the 
diagonal elements of p l 3 are similarly suppressed by at least 1/Sbh and readily 
absorbed into the normalization. 

The fact that the Hawking density matrix is diagonal can be understood by 
looking at Hawking’s derivation of the single-particle density matrix pu, as defined 
in Eq. (4.22). The relevant factor goes as ph{oj,uj') ~ /_(u;,u/) such that [1] 

oo oo 

= j dCj ( u)- 1+i = I dy = 2n6{u - J) , (5.1) 

0 —oo 

where y = In u has been used. This equation makes it clear that the matrix elements 
of p H are diagonal in frequency space. 

The off-diagonal elements of pH for different occupation numbers are also van¬ 
ishing, as follows from an estimation of the products 

/•oo /*oo 

a(J = dti otuuPu'u , atj3= dti • (5.2) 

Jo Jo 

The products af3 and aft arise in the expectation values of operators like 

(0_| ( b + ) n 6 n+2fc |0_) ~ {Pfi) n • ( a(J) k , (5.3) 

(0_| ( b + ) n+2k b n |0_) ~ [PP) n ■ (afi) k (5.4) 

7 See, however, the discussion in the Introduction about corrections arising from geometries 
without a horizon. 
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and, thus, define the elements of the density matrix between states that differ in 
total occupation number by AAf = 2k. The products a/3 and a/3 also appear in the 
computation of the generating function (A.2). Consequently, the dependence of the 
density matrix on p is only truly valid in the case when the products a/3 and a/3 are 
negligible. 

In Hawking’s Fourier-space analysis, such terms contain an integral that is similar 
to the one in Eq. (5.1) but now the two frequencies appear with the same sign, 

oo oo 

I+{uo,uo') = j dCo (£)-!+^+-') = j dy e iy(u+u') = 2 m 8(oo + oo') . (5.5) 

0 — oo 

Since both frequencies are positive, the argument of the delta function can only vanish 
if both of the frequencies vanish. This leads to one or more factors of <5 {uo + u/) when 
contracting pairs of operators in the corresponding density-matrix element, as there 
would necessarily be unequal numbers of creation and annihilation operators within 
the expectation value. 

To help reveal where corrections to off-diagonal elements can appear, let us con¬ 
sider the integration limits in Eqs. (5.1) and (5.5). These limits are an idealization, 
as a frequency can never be truly infinite nor exactly zero in physically realistic situ¬ 
ations. For instance, the upper limit on Co should be an exponentially large number, 
representing ultra-high frequencies on past null infinity that need to be red-shifted 
in order to produce the thermal radiation at future null infinity. An estimate of 
this upper limit in units of the BH mass M is ~ e T//M ~ ~ 6 s bh(o)^ 

r ~ M 3 denotes the lifetime of the BH (see, e.g., [24]). As for the lower limit on 
Co, a natural choice is to fix it as the inverse of the BH lifetime (r/M) -1 ‘S'rh^O)- 

However, in Hawking’s model, the BH is regarded as eternal, and so this limit should 
rather be the exponential of a large negative number. 

In the following, we replace the idealized limits with suitable ultraviolet and 
infrared cutoffs, e yrnax and e~ Vrnin , which are assumed to be of order y ma x ~ Umin ~ 
Sbh(0). We will further set y * = y ma x = Vmin-, as this symmetry will simplify the 
calculations without affecting the conclusions. This estimate determines the strength 
of the off-diagonal corrections of the density matrix in mode-occupation number. If, 
for some reason, the limits are such that y ma x ~ Vmin SW(0), then the original 
Hawking calculation is significantly modified and the original conclusion about the 
nature of the density matrix needs to be revised. 

The new choice of limits leads to a correction to the previous delta function, 

y* 

j dy = 2y.t^0» . (5.6) 

-y* 

The tails of this “regulated” delta function provide the off-diagonal correction of 
interest. Notice that this comes about for both I±(oo,co r ), and, hence, the off-diagonal 
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elements of the density matrix can now be non-vanishing, even when connecting 
states that differ in total occupation number. 

To better understand the implications of these off-diagonal corrections, the single¬ 
particle density matrix (and its related products) can be convolved with the same 
wave packets as in Eq. (4.21). This analysis is carried out in detail in Appendix C, 
from which there are two main conclusions: First, we find that the off-diagonal ele¬ 
ments of the single-particle density matrix p ,n -, n , = (/3/3) jrt -, n , are suppressed relative 
to the diagonal elements by a factor of p~ 2 J -1 , where /i* = y*e / (Atx 2 T) and J = j'—j. 
Second, it is found that all the elements of the products a(3, a(3 are subleading and 
suppressed as pf 2 J~p- Here, J a p = j'+j +1 and can be treated as some typical aver¬ 
age value of the frequency label j referring to the frequency range ojj E [js, ( j + l)e]. 
These products set the order of magnitude of the off-diagonal elements of the full 
density matrix. In particular, this result tells us that the leading-order contribution 
to a typical off-diagonal element of the density matrix 8 is relatively suppressed by a 
factor of iif lAMl j* 2 . Here AA f is the difference in total occupation number of 

the corresponding states and j* can be viewed as a typical value for the integer j. 

Let us elaborate on this last claim by, first, understanding the structure of the 
density matrix p. A natural way of organizing the elements of the density matrix 
which connect different Fock states is to classify them by the total number of parti¬ 
cles in each state. A density-matrix element then denotes a block of elements 
connecting states with total occupation numbers A f and A f. Thus, each element 
p\ff is a matrix of dimension djy x djg, where djg is given by the size of the Fock 
subspace spanned by the states with total occupation number equal to A f. If we 
denote the number of available frequency modes by IT, then djg coincides with the 
number of distinct elements in a symmetric tensor of rank A f in a IT-dimensional 
vector space: 

M V+ V _1 )- M 

Suppose that our interest is the strength of the density-matrix element p^j-i such 
that AA f = A f — A f is an even, non-zero integer. (There are no contributions when 
this difference is odd; cf Eqs. (5.3) and (5.4).) The strength of such an element is 
determined by a higher-order moment than the two-point function; rather, by an 
(A f + A/^-point function with A f annihilation operators and A f creation operators. 
But, because the theory of interest is non-interacting, the relative strength of such an 
element can be assessed by looking at products of two-point functions j3/3, ot(3 , a]3. 
Here, we will be able to correctly match a creation operator with an annihilation op¬ 
erator for all but of the total number of operator pairs ■ Each mismatched 

8 Typical in the sense that only an exponentially small fraction of the off-diagonal elements 
deviate from this behavior. 
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pair—two creation operators or two annihilation operators—will then contribute one 
suppression factor of pf 2 jf l . 

Formally, the density matrix can be represented in terms of the products ff, 
af, and af as 


Af = 1 -> 

tw 

0 

a/3 

0 

(af) 2 

0 . 



TV' = 2 -+ 

0 

ff 

0 

af 

0 

(af) 2 . 



Af' = 3 ->• 

af 

0 

M 

0 

af 

0 . 


(5.8) 

Af = 4^ 

0 

a/3 

0 

ff 

0 

af . 




V : 






■J 



t 

t 

t 

t 

t 

t 



N= 1, 2, 3, 4, 5, 6 

Hence, a typical density matrix element between states with a difference in the total 
occupation number A Af is suppressed by the product of factors of /i” 2 ;” 1 . 

As already stated, the off-diagonal contribution to pn —the single-particle den¬ 
sity matrix— is suppressed by a factor p~ 2 (j — j')” 1 ~ (y* £ /T)~ 2 (j — j') -1 . Since 
the final result, when traced over all the frequency modes, should not depend on the 
parameter e, we can choose it to be of the order of the BH temperature, £ ~ T. 
This amounts to saying that all the modes have the same frequency and, hence, 
j ~ 0(1). ft is then clear that the off-diagonal contribution to pu is suppressed by 
at least y~ 2 ~ 5^(0) < N ~ 2 with respect to the diagonal elements and is, therefore, 
inconsequential to the moments of pn- 

Nevertheless, the question of whether the off-diagonal elements can make a sub¬ 
stantial correction to the moments of the full Hawking density matrix pH is well 
founded, as the number of off-diagonal elements is certainly very large. To address 
this concern, let us consider the regime in which N ~ Sbh because this is when 
significant changes to the nature of the BH radiation could be expected. To estimate 
the importance of the corrections for the Renyi entropy ^ tr ^ 2 we shall replace the 
hierarchical structure of the density matrix (5.8) by a uniform estimate of the typical 
value of the off-diagonal corrections. We do so by setting |AA/"| ~ N ~ Sbh for all 
the off-diagonal elements, so that the relative suppression factor goes as (yU*J*) _T - 
The overall off-diagonal contribution to the Renyi entropy can now be calculated 
as the product of the dimensionality of the matrix d and the square of the relative 
suppression factor, d x (p 2 j*)~ N - 
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The size d of the full density matrix for N emitted particles is given by 



where d n is the size of the n-particle subspace defined in Eq. (5.7), and W is the 
number of available frequencies. Since W ~ N for BH radiation (see Subsection 4.3), 
this gives d ~ e 2N . 

Thus, the overall off-diagonal contribution, relative to that of the diagonal, is 
e 2N (yl)~ N ~ e 2N N~ 2N <C er N ~ e~ SBH , where we have again used that e ~ T and 
j* ~ 0(1), leading to the estimate njj* ~ y 2 ~ Sg H ~ N 2 . ft can be concluded 
that the off-diagonal contribution to the Hawking density matrix is exponentially 
suppressed. This conclusion agrees with the analysis of Mathur [25], who showed 
that small corrections at the level of the two-point function cannot change any of 
Hawking’s basic outcomes. The physical reason behind this result is that the thermal 
radiation at future null infinity originates from ultra-high frequency modes at past 
null infinity which are highly red-shifted due to the presence of the horizon. We 
can therefore deduce that the existence of a region of high redshift—the horizon— 
induces a characteristic exponential suppression of the off-diagonal elements in mode- 
occupation number of the full density matrix. 

Finally, let us use Eqs. (3.9) and (4.25), to determine the Renyi entropy for the 
Hawking model. As off-diagonal corrections are suppressed at least by order e~ SsH 
at all times and the diagonal is approximately uniform, the single-particle density 
matrix always has N non-zero eigenvalues, each of which is unity up to negligible 
corrections. Hence, the estimate 


(^ 2 ) Hawking — At In 3 


(5.10) 


is valid at all times. 

5.2 The Page model 

In the Page model of BH evaporation [10], the BH and the emitted radiation are 
assumed to be in a pure state in some large Hilbert space H. This space is partitioned 
into an “in” part, Hbh , representing the BH (including the ingoing radiation) and 
an “out” part, "H out , representing the outgoing Hawking radiation, so that H = 
Hbb <8> 1-Lout- The two Hilbert spaces are characterized only by their dimensionality. 

Let us label the states in Hbb by i = 1,... ,m, where m = e SBH , and the states 
in H out by A = 1,.. ., n, where n = e Srad . A state in H can then be written as | A, i). 
The working assumption is that both Hilbert spaces are large, n,m 1. Therefore, 
the Bekenstein-Hawking entropy Sbh of the BH and the radiation entropy S ra d ~ N 
(which is calculated as if the radiation was in a thermal state) characterize the 
dimensionalities of the corresponding Hilbert spaces rather than the real entropies 
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of the subsystems. Since the total state is pure, the actual entropy of the radiation 
and that of the BH are entirely due to entanglement and, thus, equal to one another 
at all times. 

The partition between “in” and “out” is meant to mimic the horizon separating 
the interior of the BH from its outside. The variation in time of the dimensions of 
the Hilbert spaces "Hbh and "H out is meant to model the evaporation of the BH, so 
that "Hbh shrinks and "H 0ut grows while the number of emitted particles increases. 
However, additional physical effects resulting from the existence of a horizon are not 
taken into account. 

The density matrix of the Page model pp is given by 


Pp = p A ,i-B,j\A,i){B,j\ . (5.11) 

The reduced density matrix for the radiation is then obtained by tracing over the 
BH Hilbert space, just as in Eq. (2.7), 

Pout = tr bh Pp • (5.12) 


Because the combined state is pure, there must be a special basis in which pp has 
only a single eigenvalue. This can be chosen, without loss of generality, to be given 
by Pi,i ; i,i- It follows that the matrix elements Pa,v,bj can be expressed in terms of 
a basis transformation matrix U acting on the original state vector 11,1). Denoting 
the resulting vector as Va, t = Ua,i- i,i, we have 


PA,i;B,j 


Ua,v,A',V PA',i'-,B',j'U^B',i'-,B,i ~ ^A,i^B,j 


(5.13) 


and 

Wp oM \B) = Y. V A,iVi,i ■ (5-14) 

i 

The vectors VA,; are meant to be treated statistically and assumed to be random 
vectors of unit size with a uniform distribution on an mu-dimensional sphere. 

The initial investigation of such a random system was conducted by Lubkin [11] 
and then improved by many subsequent investigations [26-28]. The final result is 
that the distribution of eigenvalues of p out was found to obey the Marchenko-Pastur 
(MP) law [29], which is the eigenvalue distribution for a certain ensemble of semi¬ 
positive-definite, random matrices. 

The MP distribution is generically divided into two parts, one set of large eigen¬ 
values and another set of vanishingly small eigenvalues. It depends on only two 
parameters, the total number of the eigenvalues and the participation ratio corre¬ 
sponding to the number of large eigenvalues. The distribution of the eigenvalues (up 
to a normalization convention) is given by 

Pmp( A) d\ = (max(l - c, 0) 5(A) + ^ V(A+- A)(A - A_)) d\ 

= (max(l - c, 0) <5(A) + ^-^\/4/c — (A - (1 + l/c)) 2 j dX , (5.15) 
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where 


(5.16) 


A + — (1 + 1/\/c) 2 , A_ — (1 — 1/\/c) 2 , 

and c > 0 is a parameter that determines the number of large eigenvalues. 

When c> 1, all the eigenvalues are about equal and there are no zero eigen¬ 
values whereas, for c < 1, a fraction 1 — c of the eigenvalues is vanishing. For a 
matrix of dimension n with its eigenvalues distributed according to the MP law, the 
participation ratio PR, as defined in Eq. (4.27), is given by 


PR 


nc 

1 + c ’ 


(5.17) 


The higher moments of the MP distribution are also known in terms of n and c (see 
Appendix D). 

For the Page model, in particular, 


c = 


m 

n 


(5.18) 


and so 


mn 

PR = - . 5.19 

m + n 

The early times of BH evaporation, when n <C m, correspond to c > 1 and 
PR ~ n. Essentially, all the eigenvalues are about equal and their value is determined 
by the normalization convention. This means that the Renyi entropy of the radiation 
is the same as that of the Hawking model, ( H 2 )page ~ Win3. But well after the Page 
time, when n m, one finds that c 1 and PR ~ m -C n. In this case, about 
n — m ~ n of the eigenvalues become zero and both reduced density matrices for the 
radiation and BH are becoming pure. 

The Renyi entropy of the radiation for the Page model coincides with that of 
BH and is given by 


H 2 


In PR ~ 


Inn = S rad ~ N 

\nm = S BH (N) = S BH (0)-N 


n < m , 
n > m . 


(5.20) 


This is depicted in Figure 1 as a function of the number of emitted particles. Since 
the entropy of the subsystems is entirely due to entanglement, it is symmetric if 
one exchanges the Hilbert spaces of the BH and radiation. As a consequence the 
evolution of the entropy is symmetric under reflections around the Page time, when 
the sizes of "Hbh and Rout are equal; that is, when half of the particles have been 
emitted, N Page & S BH ( 0)/2. 


5.3 The semiclassical model 

The semiclassical model [13, 14] improves upon Hawking’s framework by taking into 
account the BH’s quantum fluctuations, as well as its time dependence due to the 
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emission of the radiation. In this setup, the single-particle density matrix psc is 
no longer diagonal and, as a result, the evaporation process becomes unitary even 
though the thermal-like emission spectrum is kept. 

The basic prescription is to assign a quantum wavefunction to the collapsing 
shell of matter in Hawking’s model and then recalculate all relevant quantities as 
expectation values. The main outcome is that psc picks up off-diagonal contribu¬ 
tions that are uniform in terms of frequency but suppressed relative to the diagonal 
elements by C\I^ [ {N ) 1 where Cbr(N) = S^ ] H (N) < 1 is a classicality parameter—a 
“time-dependent h” that keeps track of how close the system is to the limit of a 
classical spacetime. 

The elements of psc do have a non-uniform suppression in terms of emission time; 
modes emitted at different times tend to decohere. Nonetheless, if the radiation is 
being regularly monitored at intervals of AN ~ \J Sbh or less, then this suppression 
can be compensated. We will specifically be considering this case, which has been 
called the “tracking case” in [17]. Hence, the off-diagonal elements of the single¬ 
particle density matrix psc can be regarded as uniform in magnitude with respect 
to both frequency and emission time. 

As explained in Subsection 4.3, p can be viewed as ai\ N x N matrix, with the 
indices running over the wave-packet modes with non-vanishing occupation number 
and with the diagonal elements given by the average occupation number for each 
mode. The elements of the semiclassical single-particle density matrix are found to 
be [13, 14] 


(. Psc)n = 1 , 

{Psc)i^j = \/CBn(N)e ldr] , (5-21) 

where the phases 0 l3 can be treated as random for most purposes (but see below). 

Let us next consider the full density matrix psc for this model. In [13] (see, 
in particular, Sect. 2.4), it was argued that the off-diagonal elements in mode- 
occupation number of the density matrix will be suppressed much in the same way 
as for the Hawking model, as discussed in Subsection 5.1. The physical reason for 
this is the same as before—the presence of the horizon implies that the frequen¬ 
cies observed at future null infinity are highly red-shifted and, thus, determine the 
magnitude of the suppression factor of the off-diagonal elements. The technical ex¬ 
planation is as follows: In the semiclassical model, the terms that lead to off-diagonal 
elements in mode-occupation number are, again, terms of the form a/3 and af3 as 
in Eq. (5.2). These contain a classical term, which is equal to the Hawking-model 
contribution, along with a semiclassical correction. Irrespective of the details, it is 
clear that the semiclassical corrections to a/3 and a(3 vanish as some power of Cbh 
in the limit Cbh —>■ 0. This is enough to guaranty the exponential suppression of the 
semiclassical contribution to the off-diagonal elements in mode-occupation number. 
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Hence, there will be a hierarchical structure in the suppression factors of the 
off-diagonal elements in mode-occupation number, as shown in (5.8). Like before, 
the off-diagonal elements can be expressed in terms of (A f + A/" / )-point functions, 
which will factorize into a product of — + 2 Jv ' / two-point functions. The strength of an 
(A/"+A/’ , )-point function will then be suppressed by a factor similar to the suppression 
factor of the Hawking model. 

In some sense, this semiclassical model can be viewed as the “middle ground” be¬ 
tween the other two models; on one hand, remaining almost thermal like Hawking’s 
with the associated hierarchical structure of the density matrix but, on the other, 
evolving over time like Page’s. At the early stages of BH evaporation, the semiclas¬ 
sical model is essentially Hawking’s plus small corrections. However, at later times 
(after the Page time), the dominant contributions will come from those elements that 
are off-diagonal in frequency. This can be attributed to the effective perturbative 
parameter being NCbu for this framework [13], as this parameter grows monoton- 
ically throughout the evaporation process and finally becomes large (.ZVC'bh > 1) 
at times later than the Page time. It is the large size of this effective perturbation 
parameter that allows the semiclassical model to evade the conclusions of Mathur 
[25]; in particular, it leads to a Renyi entropy for the radiation that is decreasing in 
time (see below). 


5.3.1 Semiclassical Renyi entropy 

The objective here is to estimate the Renyi entropy of the BH radiation H 2 (N) using 
the single-particle density matrix psc as specified in Eq. (5.21). However, we are 
currently lacking a precise quantitative knowledge of p$c- bi particular, in order to 
use Eq. (3.9) for H 2 (N), one needs to know the eigenvalue distribution P(X) of psc'- 


H 2 (N) = (tr [In (1 + 2p)])p ( jq — N In (1 + 2X)P(X)d\ . (5.22) 

Jo 

This deficiency can be attributed to our lack of specification as to the exact nature 
of the phases e l9ij . These phases are not entirely random but rather dictated by the 
effect of the quantum horizon on correlations in the emitted radiation [14]. In this 
subsection, we will instead call upon the known properties of psc hi order to draw 
conclusions about the entropy of the radiation in the semiclassical model. 

We expect that psc has some number of large eigenvalues with the rest vanishing, 
whereby the large eigenvalues are distributed about the mean in a way that depends 
on their total number. This is most clearly seen from the participation ratio PR of 
the single-particle density matrix psc■ The participation ratio provides, for a given 
matrix, a measure of the number of non-vanishing eigenvalues (e.g., [23]). For psc, 
it is given by [17] 


(tr Psc) 2 _ N _ N(S bh (0)-N) 
tr p 2 S c 1 + IVCbh S bh { 0) 


(5.23) 
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Then, since tr psc = N (as shown in Subsection 4.3), the average value of a non¬ 
vanishing eigenvalue can be estimated as 

J=tI Ff = TR = 1 + NC ™' < 6 ' 24) 

This finding reveals that, up to the Page time (IVCbh < 1), there are about 
N eigenvalues of order one. The properties of the emitted radiation in this regime 
coincide with the predictions of the Hawking model: Each of the N available modes 
has an occupation number of order one and the radiation entropy grows linearly with 
the number of emitted particles. 

On the other hand, after the Page time (IVCbh > 1), the participation ra¬ 
tio (5.23) stops growing and eventually becomes much smaller than N, so that the 
average A is much larger than one. This signals that the matrix p$c has a small 
fraction of large eigenvalues, with the majority of eigenvalues vanishing. Physically, 
it means that, out of the N available radiation modes, there is only a small fraction 
with large occupation numbers. In this regime, the entropy of the radiation (5.22) 
is dominated by the contribution from the occupied modes and is smaller than the 
entropy of the thermal radiation (5.10) which is predicted by the original Hawking 
model: 


H 2 (N) ~ PR In (1 + 2A) ~ PR In 



<C N , 


for PR <C N . 


(5.25) 


Therefore, the known properties of the single-particle density matrix psc guarantee 
that the radiation entropy in the semiclassical model does not grow steadily with 
number of emitted particles, as predicted by the Hawking computation, but rather 
starts to decrease after the Page time. 

In order to illustrate the evolution of the radiation entropy in the semiclassi¬ 
cal model, let us consider two examples for the eigenvalue distribution of p S c, both 
having the property that the participation ratio can be much smaller than the dimen¬ 
sionality of the matrix. First, we consider the simplest case when psc has exactly PR 
non-vanishing eigenvalues of magnitude A and the remaining (N — PR) eigenvalues 
are exactly vanishing. The corresponding distribution function reads 


P{ A) 


N — PR PR 

<5(A) + — 5 


N 


N 



(5.26) 


where the participation ratio is that of Eq. (5.23). The Renyi entropy in such a case 
is given by 


H 2 (N) = PR In (1 + 2A) = 
= Sbh( 0) 


N 


1 + IVCbh 
N 


Sbh( 0) 


In (3 + 2 IVCbh 
N 


Sbh( 0) 


In 1 + 


1-N/S B h(0) 


. (5.27) 
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N/S bh (0) 


Figure 1 . Dependence of the Renyi entropy H 2 of the BH radiation on the number of 
emitted particles N for the semiclassical model in Eq.(5.21) for two eigenvalue distributions. 
First, assuming that the non-vanishing eigenvalues of the two-point function are distributed 
according to Eq. (5.26) (solid, thick orange) and then to the Marchenko-Pastur law in 
Eq. (5.15) (solid, thin green). For comparison, the dashed blue line depicts the prediction 
of the Page model from Eq. (5.20). 


The evolution of this radiation entropy with the number of emitted particles is pre¬ 
sented in Figure 1, which shows the predicted decline after the Page time. 

Another eigenvalue distribution that illustrates the properties of psc is the MP 
distribution (5.15), as already discussed in the context of the Page model. The 
crucial difference from the Page model itself is that we now use the MP distribution 
to model the single-particle density matrix p S c and not the reduced density matrix 
which represents the state of the emitted radiation. Requiring the participation ratio 
of the MP-distributed matrix (5.17) to be equal to the semiclassical prediction (5.23), 
one finds that the parameter c of the MP distribution is given by 


1 

IYCbh 


(5.28) 


The Renyi entropy of the radiated particles can then be obtained by plugging the 
MP distribution (5.15) into Eq. (5.22). The evolution of the entropy is shown in 
Figure 1. One can observe that the entropy is decreasing after the Page time, just 
as expected for the semiclassical model. 

Although these two examples were used for illustrative purposes, we do expect 
the eigenvalue distribution for psc to be similar in some aspects to the MP dis¬ 
tribution. The moments of the distribution function of the eigenvalues of psc are 
calculated and compared to the moments of the MP distribution in Appendix D. In 
spite of differing in details, the two sets of moments follow a similar pattern. 
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Having established that the Renyi entropy of the radiation in the semiclassical 
model starts to decrease after the Page time, let us speculate on a viable form for its 
functional dependence on the number of emitted particles. What is known about the 
semiclassical model is that H 2 (psc ) (1) must tend to the Hawking value of Wind, as 
in Eq. (5.10), at early times and (2) must be symmetric under an exchange between 
the radiation and the BH, since the two reduced density matrices necessarily share 
the same set of eigenvalues. In particular, all physical quantities are required to 
be symmetric under the exchange of N O Sbh(N), as long as these quantities 
characterize the sizes of their respective Hilbert spaces "H out and 'Hbh- 

The participation ratio (5.23) is manifestly symmetric in just this way, 


PR(psc) 


N 

1 + NC bh (N) 


NS bh (N ) 
Sbh(N) + N 


(5.29) 


It is also of the same form as the participation ratio of the MP distribution (5.17), 
given that one identifies the parameters n and m of the Page model to be N and 
Sbh(N), respectively. 9 This may seem peculiar because, as discussed in Appendix D, 
the actual eigenvalue distribution of p$c differs from the MP law. But, as far as we 
can tell, simplest non-trivial, symmetric function which tends to N 

for N <C Sbh (0), suggesting that such an identification is valid. 



n/s bh (0) 


Figure 2. Dependence of the Renyi entropy H 2 of the BH radiation on the number of 
emitted particles N. The solid purple line depicts the prediction of the semiclassical model 
from Eq. (5.30) and the dashed blue line depicts the prediction of the Page model from 
Eq. (5.20). 


9 Note that, in the Page model, n and m are the dimensions of the reduced density matrices p ou t 
and pp, and are exponentially large numbers, n ~ e N and m ~ s s bh(n) _ 
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Motivated by these observations, our expectation is that the semiclassical Renyi 
entropy can be expressed in terms of the participation ratio as follows: 


H'lipsc) 


PRiPsc) 


In 3 + c_ i 


1 

PRiPsc) 


1 

+ c - 2 ( pfl ( Psc )) 2 + "'. 


(5.30) 


where the coefficients are dimensionless numbers that are determined by the higher 
moments of the eigenvalue distribution. To see that Co = In 3 is indeed the correct 
leading-order coefficient, consider that, at early times, N ~ PR to leading order in 
a 1/PR expansion. Thus, the average eigenvalue A = ~ 1 and the logarithm in 

Hi then expands as In (1 + 2A) ~ In 3. 

The dependence (5.30) of the Renyi entropy on the number of emitted particles 
for a semiclassical BH is presented in Figure 2. Because the Renyi entropy depends 
only on the participation ratio 1+j ^; BH , one can determine the evolution of the for¬ 
mer from that of the latter. The participation ratio steadily grows until reaching 
a maximum at the Page time and then steadily declines for the remainder of the 
evaporation process. This is qualitatively similar behavior to the Renyi entropy for 
the Page model although, in the semiclassical model, the rate of purification—which 
is the rate of deviation from the linearly growing result of the Hawking’s model—is 
actually faster than the Page-model rate [17]. 


6 Conclusion 

We have exposed the relationship between the density matrix of the outgoing radia¬ 
tion from a BH and the corresponding single-particle density matrix in the case that 
the density matrix is approximately diagonal in mode-occupation number. It was 
then shown that the presence of the horizon leads to a high suppression of the off- 
diagonal elements in mode-occupation number of the density matrix for the emitted 
radiation. We have therefore concluded that the density matrix is approximately 
diagonal in mode-occupation number at all times. We have also shown how to reg¬ 
ularize the infinities which arise as a consequence of Hawking’s idealized picture of 
an eternal BH. It was explained how the state of the emitted radiation from a BH is 
different from that of a standard black body of the same temperature. This analysis 
was then applied to three models of BH evaporation as a means of clarifying their 
differences and contrasting their main features. Let us briefly summarize the main 
observations: 

The Hawking model: 

The off-diagonal elements of the single-particle density matrix are highly sup¬ 
pressed for different frequencies and emission times, and those of the full density 
matrix are exponentially suppressed for differences in occupation numbers. The 
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Renyi entropy (4.20) scales with the number of emitted particles N at all times, and 
so there is no possibility for the radiation to be purified. 

The Page Model: 

For a randomly chosen basis whose relationship to Hawking’s basis is unspeci¬ 
fied, the off-diagonal elements of the density matrix are uniform in magnitude with 
random phases. There is no hierarchy between different elements of the matrix. The 
Renyi entropy of the radiation (5.20) increases linearly as N until the Page time and 
decreases linearly as Sbh( 0) — N afterwards. Hence, purification is inevitable. 

The semiclassical model: 

The off-diagonal elements of the single-particle density matrix are suppressed 
by a power of S^J 2 in a uniform way, whereas those of the density matrix are ex¬ 
ponentially suppressed for different mode-occupation numbers in a similar way to 
the Hawking model. At the Page time, the contribution from the off-diagonal el¬ 
ements with regard to frequency and emission time becomes significant because of 
a perturbative parameter that grows monotonically with time. The Renyi entropy 
(5.30) scales with the participation ratio of the single-particle density matrix, which 
increases until the Page time and decreases thereafter. This suggests that the radia¬ 
tion starts to purify and does so at a rate which is faster than that of the Page model. 
The question of how unitarity is recovered for a process of gravitational collapse in 
the semiclassical model is discussed in detail in a companion paper [17]. 
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A Expression for the full density matrix in terms of //• 

In order to find an expression for the density matrix, it is useful to calculate a vacuum 
expectation value; namely, 

(0_|e 6 ^e A ^|0_) . (A.l) 

The resulting function of /P and A j can be used as a generating function for the 
expectation values of any normal-ordered powers of the creation and annihilation 


operators bf and IP. The expectation value (A.l) can be calculated by using (3.2) 
and applying the Baker-Campbell-Hausdorf (BCH) formula for the exponents: 


(Q^\e b i ^e Xjbi |0_) = e~^ il3ijaj e~ Xj0jka ^ e Aia A afc |o_) e~^ ia * ^ kjpk e~^ Xiaik ^ jkXj 

(A.2) 


The exponential factors on the outside of the average arise from the commutators. 
By using Hawking’s expressions in [2] for the Bogolyubov coefficients a t j and fy, one 
can show that 

a/Ay = 0 , a\F k = 0 , (A.3) 

and so the external exponentials are equal to unity as their exponents are vanishing. 

Relations (A.3) also ensure that only the matrix elements of p out between states 
with the same total occupation number are non-vanishing. But, as discussed in 
Section 5, these relations hold only approximately for realistic BHs. 

Returning to Eq. (A.2), one can see that the two outer factors within the expec¬ 
tation value disappear after acting on the vacuum |0_). The remaining expression 
can be further simplified by applying the BCH relation once again: 

^0 | e -/hft^ e -A^' fc a+|0 \ = g| e -A^' fc a+e—\ _ gAi/3 lfc /3jfc/h 

(A.4) 

Finally, the answer for the generating function (A.l) reads 

(O-le^V^lO.) = e Xip > j , (A.5) 


where p* is the single-particle density matrix as defined in Eq. (3.3). For Hawking, 
p l - is a diagonal matrix but, in more general setups, it can also have non-zero off- 
diagonal elements. In any case, p x 3 is Hermitian and can thus be diagonalized by a 
unitary transformation. 

We next look for a closed expression for the density matrix p° ut in terms of the 
6’s. In analogy with the harmonic oscillator, let us try the ansatz 


P 


out 


1 -fe+fn.y 

Z 


(A.6) 


where H is some c-number Hermitian matrix and Z is the normalization factor. The 
matrix H can be diagonalized by a unitary transformation, O = U + DU, with D 
being a real diagonal matrix. This allows us to introduce new annihilation operators 
d 7 = U\b l for which the quadratic form b + flb = b + U + DUb = d + Dd is diagonal. 
Such a transformation is canonical, so that {df, di} are also a good set of creation 
and annihilation operators. 

Moreover, since this transformation does not mix the creation and annihilation 
operators between each other, the subspaces with a fixed total occupancy of b- and 
d-particles coincide, and the traces over the Fock spaces of b- and d-particles are 
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equal. In particular, one can find the normalization Z by summing over the states 
with a fixed occupation number of rf-particles |rq): 


Z — tr 




= tr 


a- d t Di i d ' 


n 


1 — e 


-D\ 


iie 

i m=0 


= det 


nAe 




1 — e 


—D 


= det 


1 — e 


-n 


• (A.7) 


The last equality is due to the invariance of the determinant under unitary transfor¬ 
mations of the matrix. 

According to ansatz (3.4), the corresponding generating function (A.l) is 



e b+ ^e Xb e~ b+nb 


(A.8) 


If it is possible to choose so that the result matches (A.5), then ansatz (A.6) gives 
the correct expression for the density matrix. 

As before, the generating function (A.8) can be calculated by tracing in the 
d-particle basis: 


1 

Z 


tr 

' e b+n e \b e -b+Ub 

= -tl \e d+U ^e X U + d e ~d+D d 1 



z L J 


1 oo 

rr 1. y^{ ni \e d n u ^ e ^ u+ )^ e -4D^\ ni) ? (A 9) 


rii=0 


where each of the factors in the bottom line corresponds to a single oscillator. 
In order to evaluate the above expression, we will use the identity 


(i 


—D 


) tr e «V e Aa e -Da+a 


= exp 


A n 


oD 


(A.10) 


which is derived in [30]. (For a web resource, see, e.g., [31].) The expectation 
value (A.8) then reads 


(e b+ »e Xb h = 


)p = II Gxp y 


/ (A U+MUfi? 


e U i — 1 


= exp ( A U 


exp ( ^ 
1 


(A U+Mutf} _ 


D l i 

e u i — 1 


e D - 1 


U n ) — exp I A 


e Q — 1 


h) • (A.11) 


Comparing the last expression to (A. 5), we find that the right-hand sides coincide 
if the matrix is related to p as 


1 


P = 


or e x = 1 + 


P 


(A. 12) 
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B Grey-body factors 


From the properties of the grey-body factors, one can see that the state of the field 
outside a large BH is quite different from the usual thermal state in an empty box, 


p = T?e ,3H . The actual state of the radiation is that of a field in thermal equilibrium 


with the BH, which is very selective in what it absorbs and emits. For example, the 


occupation number of the modes with high angular momentum will be much less 


than what one would naively expect from just the Boltzmann suppression. 

One way to present this modification away from thermality is to notice that, by 
definition, the grey-body factor T u i m is equal to the modulus square of the scattering 
amplitude of the corresponding mode on the BH. Therefore, it is physically intuitive 
to introduce a corresponding effective scattering cross-section (see, e.g., Ref. [22]), 


(B.l) 



The rate of emission can then be written as 


A"a jlm f duj Im 

5t L 2n 2 e u / T — 1 


f 

Jo 


(B.2) 


This is a three-dimensional thermal emission rate, where the geometric area of the 


emitter is different for each mode and is given by the absorption cross-section. In this 


language, the suppression of the high-/ modes comes from the fact that the effective 
area of the BH as seen by them is much smaller than its geometric area. 

Another way to see the difference from the usual thermal state is to notice that 


the density matrix of the equilibrium state has the form p = Z l e b+nb with fl = c o/T 


only if T uim = 1. Indeed, the relation (3.5) between the single-particle density matrix 


p = and the matrix H implies that 



(B.3) 


Therefore, the density matrix for the radiation takes the form [cf. Eq. (2.10)] 



(B-4) 


For each given mode, the ratio of the probability of having (N + 1) particles to 
the probability of having N particles is constant and given by 


p(n +1) _ r uim 


(B.5) 


P(N) e“l T - (1 - r„ im ) ' 


Meaning that, apart from the overall Boltzmann suppression factor, the emitted 
particle has a probability of (1 — T) to be scattered back into the BH, 


P(N + 1) = e- u/T [T P(N) + (1 - r) P(N + 1)] . 


(B.6) 
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C Off-diagonal corrections to the Hawking model 

In order to estimate the off-diagonal corrections to the Hawking density matrix, we 
first need to consider the corrections to the products 

P ik P jk , a/Pkj , a{P jk , (C.l) 

where the second index is summed over. The first product defines the standard 
single-particle density matrix p* = f3 lk /3jk and, in the case of Hawking’s calculation, 
is diagonal and completely determines the density matrix of the outgoing radiation 
(3.4). The products a/3 and af3 arise in the expectation values of operators like 

(0_| {b + ) n b n+2k |0_) ~ {Pfi) n ■ (a(3) k , (C.2) 

(0_| [b + ) n+2k 6 n |0_) ~ (£0)" • (a/3) k (C.3) 

and, thus, define the elements of the density matrix between states that differ in 
total occupation number, A J\f = 2k. 

The products a/3 and a/3 also appear during the computation of the generating 
function (A.2). Consequently, the sole dependence of the density matrix on the 
single-particle density matrix p* is only truly valid in the case when the products a/3 
and a/3 vanish. This assumption is true for Hawking’s idealized calculation but, as 
shown below, modified for physically realistic BHs. 

In order to find the corrections to the products (C.l), we rewrite the definition 
of the coefficients aw and /3 UUJ i in the basis of Fourier modes, 

= if a (u>) (cu)- 1/2+ ^ e^-^o , (C.4) 

/U = -ifpiu) (A)" 1/2+ ^ e^ + ^vo } (C.5) 


where 

fa M = (l - —\ , f p (u) = f a (u)e~ , (C.6) 

Z7T y/CJ y K J 

and Vo denotes the position of the BH horizon in advanced time, k = 2ttT is the BH 
surface gravity and t u is the transmission coefficient for which \tj\~ = T w . 

The products (C.l) can be rewritten in the wave-packet basis of Section 4. For 
instance, 

r(J+ i)e rti'+ 1) £ 

(ppy n =e-' due 2 "™/’ (C.7) 

Jje J j'e 

where 

/»oo 

P0p(u,<jj')= / duPuuPw'Cj- (C.8) 

Jo 
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There are similar expressions for (a:/?) ■ ■/„/ and ( af3y nj " . The product (C. 8 ) and 
its analogues can be expressed in terms of the Fourier-mode basis as 


Ppp(u,U}') 


— U))vo 

r+oo 

1 dye iyl ^'-“\ 

(C.9) 





r+oo 


p a p{u,u') 

= fa 1 

+)M+)e^ 

'+u)vq 

/ dye iy/K{uj ' +u) , 

(C.10) 





r+oo 


Papfau') 

= /+ 

+)M+)e~’ ( 

Ul'+u)vQ 

/ dy e -+/«+'+^) , 

J — OO 

(C.ll) 


where the integration variable has been changed to y — ln(av). 

Substituting the previous set of relations into the corresponding expressions for 
the products in the wave-packet basis and then integrating over the frequencies oj 
and u /, we arrive at 


m‘ n iv 

( a fi)j n j'n' 

/-n\jnj'n' 

\<*P) 


2k f + °° g27T in(j'-j) 

fp (^') fpfi) — dfi -—-rw sill 2 (vr/i) , 

K J -oo (p — n)(/i — n') 

2 k /' +0 ° p^mO'+j+i) 

fa (+') f p {Q)— / dfi - - — sin 2 (vr/i) , 

2k r +oa p—27ri/i(j , +j-|-l) 

fa (&') ff 3 (w) — / dfi -—-— sin 2 (vr/i) , 

7T J -oo (/i-n)(/i-n') 


(C.12) 

(C.13) 

(C.14) 


where a) = £ (j + 1 / 2 ) and u’ = e {j’ + 1 / 2 ) are the mean values of the frequencies 
in the given range of integration, n = n — no, n' = n' — no with no defined by 
no = 2vrno£ _1 , and the integration variable has been changed to n — ye/{2n k). We 
see that the three integrals differ only by the power in the exponent, and so it is 
convenient to define 


Jpp = f — j 1 Ja/3 = j' + j + 1 ) Ja/3 = ~{f + 3 + 1) ■ (C.15) 


Let us now introduce infrared and ultraviolet cutoffs for the integration range of 
the frequencies, as discussed in the Subsection 5.1, 



-+ 



(C.16) 


Due to the similarity of the three integrals above, we need only to consider one, 


/ +M* g27T IflJ 

(p. - n)(„ - n') Si “+ ri ’ (C ' 17) 

and then interpret the results for the different choices of J in Eq. (C.15). The integral 
does need, however, to be treated differently for the cases with An = n' — n = 0 and 
An 7 ^ 0, as well as for the cases with J = 0 and J / 0. 
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C.l Off-diagonal elements in frequency 

When An = 0, the integral (C.17) becomes 


i°(j) = 




d/j, ■ 


0 2nifiJ 


■ sin 2 (7r/i) . 


j —ju* 

For the case J = 0, the previous integral can be evaluated to give 


(C.18) 


/o(0) = 


sin 2 ( 717 ^) sin 2 ( 7171 *) 


+ 7rSi(27r(/z* — ft)) + 7rSi(27r(/i* — ft)) , (C.19) 


/i* - ft /i* T ft 

where Si(a;) = dt sin t/t is the sine integral. It has the large-£ expansion 

, 7 r cosx sin x 

Si (x) = - -- + °( x ) • 

A X X A 

Hence, in the limit when /i* ^ n, the integral becomes 

T _ „2 1 sin(27r/r*) //Vi _ 3N 

Jo(0) — 7T--- 2 -S 0(n* ) . 

h* 2 tt/^ 


When J 4 0, the integral can be split as 


I 0 (J) = x h(J) - -h(J - 1 ) - -h(J + 1 ) 


where 


h(J) = 




0 2lTl[lJ 


'-,u - n) 2 ' 


(C.20) 

(C.21) 

(C.22) 

(C.23) 


We hnd that the last integral is expressible as 


g27r i/itJ g—27r ifitJ 

h(J) = ---—- AniJffi (— AiriJ (/ 1 * — ft)) + 27tzJE 1 (27uJ(/i* + ft)) , 


(C.24) 


[i* — n /i* + n 

where E\(z) = dte~ zt /t is the exponential integral with the large-x expansion 

(C.25) 


Ei(zx) = e-“ ( - + \ + 0(x-*] 


IX X* 


I 11 the limit when /i* ft, the function h(J) then becomes 


-£ + 0 (/i 




J = 0 , 


h(J) = 


+ J* o. 


(C.26) 


fit 2-ITfltJ 

The integral / 0 (J) at a given value of J can then be evaluated by combining 
Eqs. (C.21), (C.22) and (C.26). For example, J 0 (l) works out to be 


J 0 (l) = - h{ 1 ) - -h( 0 ) - -h{ 2 ) 


1 f sin(27T/i*) 1 1 sin(47T/i* 


2 V K /4 


/i* 




(C.27) 
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C.2 Off-diagonal elements in mode-occupation number 

When An X 0, the integral (C.17) can be split as 


Ia„(J) = l (g(J) - \g(J - i) - \g(J +1) ] , 


(C.28) 


where 


9(J) = 


1 

An 




djji 


D 27riJ 


— (n 


(C.29) 


--/i \V-n' 

If J = 0, the integral can be computed exactly. For J ^ 0, it is possible to 
redefine the integration variables in the two summands so as to move the dependence 
on the expansion parameter n* from the limits to the integrand: 


g(J) = 


2cos(2vr/i*J) /'" ^ e~ 2nitJ | 2i sin(2vr/i, J) f h ^ e~ zmtJ t 

t \ 2 n* 


Ann* 


dt ■ 


1 - ( - 


+ ■ 


An/j, 


— 27 ritJ 

dt —-— : —— — . (C.30) 


1 - - 


One can then expand the integrand in the limit when /j,* n, n! and evaluate the 
integral to obtain 


g(J) = 


—~jX + 3 )) J — 0, 

+ ofc 3 ), J* o. 


(C.31) 


Comparing g( J) to the function h(J) in Eq. (C.26), we see that the two functions 
are equivalent, g(J) = h(J). Hence, for J X 0, the two cases are coincident, 


Iau(J) = W) , <7^0. 


(C.32) 


On the other hand, at J = 0, the integral I/± n differs from 7 q( 0) in Eq. (C.21) 
only by a zeroth-order, diagonal term: 

W0) = 5 (h( 0) - \k- 1) - ift(l) 


1 sin(27r/i* 
n* 2vr/i2 
= J 0 (0) - 7T 2 . 


o(/i 


- 3 \ 


(C.33) 


Hence, the corrections in both cases, An = 0 and An X 0) coincide and are 
uniform in An. 


C.3 Final Result 

Up to the leading order in fi* 1 , the integral (C.17) then becomes 

iAn(J) — ^dAnfibjfl + — <5/,o + -<5/, 1 + -Sj ,+ 0(n* 2 ) j (C.34) 
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where the second-order corrections are subleading for J = 0, ±1. In the other cases— 
he., when J 7 ^ 0, ±1— it is the corrections at second order in /i ^ 1 that are dominant: 


Iau(J) 


1 / sin(27r/i* J) 
/i* \ 27r/i*J 


l sin(27r^(J-1)) l sin(27r^(J +1)) \ _ 3 

2 2tt//*(J-1) 2 2tt/ 7*(J + 1) J + l/ "* 


(C.35) 

In order to interpret these results as corrections to the density-matrix elements, 
one has to employ the different definitions of J in Eq. (C.15) for the different cases. 


Interpretation 

The matrix elements of the product (/3/3)‘ ?n , n , are proportional to 


W" ,, n , - I An {J m ) , J m = j' - J . (C.36) 


Since the corrections are found to be uniform in An = r! — n, we can represent the 
matrix (j3/3y n , n , as a tensor product of j’s and rC s. Then each matrix element that 
is labeled by j and j 1 is itself a uniform matrix in terms of n and n'. I 11 this way, the 
product can be represented in the following form: 


m 


J n 


3 n' 


/ 1 

1 



1 

1 

1 

t: I 

* 1 

1 


ifij 

1 


J_ 

/n 

1 A. 

11 * 




(C.37) 


where each depicted entry represents a uniform block and only the order of magnitude 
of the leading term is shown. 

For the products a/3 and a(3, the corresponding index is J a p = ±(j' + j + 1), 
respectively. Since both j and j' are positive, these products have no zeroth-order 
contributions and only the matrix elements with j = j' = 0 have corrections of 
the order 0( l//i*). All the other entries of the matrices a/3 and a[3 receive only 
second-order corrections O (l/(/i^ J)). 

How these findings impact upon the off-diagonal elements of the full density 
matrix is discussed in Subsection 5.1 of the main text. 


D Higher moments of the Marchenko—Pastur and psc dis¬ 
tributions 

Here, we will quantify more precisely the differences between the Marchenko-Pastur 
(MP) distribution and the eigenvalue distribution of psc■ This entails assigning the 


36 









distributions the same participation ratio and then determining how their higher 
moments are different. 

Comparing their respective participation ratios (with n = N) in Eqs. (5.17) 
and (5.23), one can readily identify the parameter c for the semiclassical two-point 
function, 

" = NCbs(N) ' (EU) 

Let us consider the case NCbh > 1- This is really the regime of interest, since 
smaller values of IVCbh correspond to the case in which the Hawking model is valid 
(up to small corrections). We will further assume that the phases of the off-diagonal 
terms in psc can be treated as random. One can then estimate the higher moments 
of psc up to combinatorial factors and sub-leading terms in small 1/NCbk , which 
leads to 


tr (/4c) 

(tr psc) 2p 

tr {psc ' 1 ) 

(tr psc) 2p+1 


N? +l (Cbh) p (Cbh) p 
7V 2 p Np- 1 ’ 

jV p+1 (Cbh) p N (C BH ) p 
N 2 p+ i Np - 1 


(D.2) 

(D.3) 


The basic idea behind these estimates is that the off-diagonal parts of the matri¬ 
ces are dominant when IVCbh ^ 1 and the randomness of the phases requires these 
parts to sum coherently (he., restricted to sums of the form JA. MijMji). A simple 
example should suffice to illustrate the point. Let 7 and 77 be the diagonal and off- 
diagonal parts respectively of a matrix p. Applying the rule of coherent summation 
and, otherwise, insisting on the maximum power of 7 , we have, for the p — 2 case, 


trp 4 = y p l3 p jk p kt pa ~ y raj Vp rju pu , (D.4) 

ijkl ijl 

^ P ^ Pij Pjk Pkl Plm Pmi ^ ^ ^ Pij Vji Pil Pli ”T ii • (D.5) 

ijklm ijl 


As one can see, each of these traces results in 3 = p + 1 independent summations, 
a trend which continues for any value of p. This accounts for the factors of A ^ p+1 
in Eqs. (D.2) and (D.3); the rest is determined by the magnitude of the elements. 
These results indicate that the higher moments of the eigenvalue distribution are 
determined by an expansion in \JC-qu- 

The moments of the MP distribution, on the other hand, are expressed as a 
power series in c and not square roots thereof. Indeed, for the same conventions and 
the same regime of small c = 1/(A^Cbh), the MP distribution would yield for the 
high moments 


(tr p M p) n \cj \N) 


(Cbh) 


n— 1 


for n>l. 


(D. 6 ) 
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This follows from the observation that the MP distribution (5.15) has, for small values 
of c, about Nc ~ 1/Cbh large and (roughly) equal-valued eigenvalues A ~ 1/c. 

This discrepancy between the MP distribution and the semiclassical distribution 
is a consequence of the square root of C B h appearing in the off-diagonal elements 
of psc■ Hence, psc does not precisely conform to an MP distribution nor should 
it necessarily be expected to. However, when 1VC B h ~ 1, both expressions for the 
moments of the distributions scale in the same way, (Cbh) ti 1 - Therefore, we 

do expect that the two distributions share the same general features; in particular, 
once c = 1 (IVCsh = 1 ) is reached, both eigenvalue distributions begin to develop 
support for zero eigenvalues. 

Yet, when NCb h 1 is true, the higher moments of the eigenvalue distribution 
of psc are much more suppressed than those of the MP-distribution. For instance, 
the n th semiclassical moment is smaller by a relative factor 1/ (NCbh) 71 ^ 2 than its 
MP counterpart. Nevertheless, we do expect that, in this case, both distributions 
have only a few large eigenvalues but apparently differ in the detail. It would be 
interesting to find out what is the actual eigenvalue distribution of the semiclassical 
matrix. 
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